%% task6_2_2_3.m

clear;
load("../data/Guitar.mat", "wave2proc");

% plot(wave2proc);
Samp = 8000;        % 采样率

% 使用傅里叶级数近似分析
OnePeriod = wave2proc(1:25);        % 近似认为25个数据点为一周期，周期T=(1/8000) * 25
FourSe = FourierSeries(OnePeriod, 25);
bF = Samp / 25;         % 基频
fprintf("近似处理: \n基频为: %f\n", bF);
% 最近的音调为e
for i = 1:5
    fprintf("谐波分量%d: 幅值: %f\n          相位: %fpi\n", i, abs(FourSe(i+1)), angle(FourSe(i+1))/pi);
end

% 对整个信号求傅里叶变换
L = length(wave2proc);
Y = fft(wave2proc);
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Samp*(0:(L/2))/L;
figure
plot(f,P1) 
title("wave2proc的频谱分析")
xlabel("f (Hz)")
ylabel("|P1(f)|")
% 观察图像，对应基频为e，C大调的Mi

% 时域重复10次，观察频域变化
wave2long = repmat(wave2proc, 10);
Llong = length(wave2long);
Ylong = fft(wave2long);
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 Llong 计算单侧频谱 P1。
P2long = abs(Ylong/Llong);
P1long = P2long(1:Llong/2+1);
P1long(2:end-1) = 2*P1long(2:end-1);

f = Samp*(0:(Llong/2))/Llong;
figure
plot(f,P1long) 
title("wave2long的频谱分析")
xlabel("f (Hz)")
ylabel("|P1long(f)|")

[pks, locs] = enve(P1long, 3);
hold on
plot((locs-1)*Samp/Llong, pks);
% xline(locs*Samp/Llong);

overtone = [];
lo = [];
j = 1;
for i = 1 : length(pks)-1
    if abs(pks(i+1)-pks(i)) > 0.0025 && pks(i+1) > 0.001
        overtone(j) = pks(i+1);
        lo(j) = locs(i+1);
        j = j+1;
    end
end
xline((lo-1)*Samp/Llong);
plot((lo-1)*Samp/Llong, overtone, 'b');
s = sum(overtone);
overtone = overtone / s;        % 归一化
funFre = (lo(1)-1)*Samp/Llong;  % fundamentalFrequency

save("..\data\wave2fft.mat", "overtone", "funFre");

% 迭代计算包络
function [pks, locs] = enve(data, iter)
    if iter > 1
        [pks, l1] = findpeaks(data);
        [pks, l2] = enve(pks, iter-1);
        locs = l1(l2);
    else
        [pks, locs] = findpeaks(data);
    end
end

function fres = FourierSeries(data_se, n_max)
% 傅里叶级数
% data一周期内离散数据集，T 周期，Ts抽样时间间隔，n_max 级数展开至第几项
% fres：1 x （n_max + 1）复数数组。

    fres = complex(zeros([1 n_max]));
    for k = 0 : n_max-1
        for n = 0 : length(data_se)-1
            fres(k+1) = fres(k+1) + data_se(n+1) * exp(-n*k*2*pi/length(data_se)*1j);
        end
    end

end